%%%%%%%%%%%%
% Summary of pricing and return results
% Liuren Wu, liurenwu@gmail.com
%%%%%%%%%%%%
format compact;    clearvars;
load('./data/Momentsw.mat');
load('./data/PricingResultsFull.mat'); %full sample results

paperfd='./temp/'; %save temporay results
portfolioweight=1;
marketpricing=1;
pricingperformance=1;
returnfactormodel=1;

pdeltap=pdelta*100;
matm=mat*12;
xm=repmat(pdeltap(:)',nm,1);
ym=repmat(matm(:),1,nk);


xpname={'\eta_V','\eta_G','\eta_O','\eta_A'};
pname={'Vega','Gamma','Volga','Vanna','Error'};
zname={'Vega','Gamma','Volga','Vanna'};nz=4;
yeara=num2str([1996:3:2023]');
if portfolioweight
    mwx=nanmean(wx,3);
    ppd=pdelta*100;
    for c=1:2
        for k=1:np
            w=reshape(mwx(:,k,1,c),nm,nk)';
            figure(k)
            plot(ppd, w(:,1),ppd, w(:,2),'--',ppd, w(:,3),'-.',ppd, w(:,4),':o',ppd, w(:,5),':','LineWidth',4,'MarkerSize',20)
            xlabel('Put delta, %')
            ylabel('Portfolio weight')
            set(gca,'Box','on','LineWidth',4,'FontSize', 25)
            print('-depsc','-r70', [paperfd,'fig',pname{k},'xweight',currency{c},'.eps'])
        end

    end
end
if marketpricing
    disp('Market pricing estimates');
    x=[RRP(:,:,1), RRP(:,:,2)];
    Tn=sum(isfinite(x));
    a=[nanmean(x),
        nanstd(x),
        sqrt(Tn).*nanmean(x)./neweystd(x),
        ];
    corx=100*diag(corr(RRP(:,:,1), RRP(:,:,2),'rows','pairwise'))';

    fprintf(1,' \n')
    for j=1:np
        fprintf(1, ['%20s '], pname{j});
        fprintf(1, [' &  %8.2f  &  %8.2f  &  %8.2f       &&  %8.2f  &  %8.2f  &  %8.2f   &&  %8.2f      \\\\ \n'], ...
            [a(:,j)',a(:,np+j)',corx(j)]')
    end
    mineta=[-20,-2,-2,-2];
    maxeta=[40,6,10,2];
    for j=1:np
        eta=[RRP(:,j,1), RRP(:,j,2)];
        figure(j);clf;
        plot(wdate, eta(:,1),wdate,eta(:,2),'-.',wdate,wdate*0,':','LineWidth',4)
        xlabel('Year','FontSize',20);
        ylabel(xpname{j},'FontSize',20);
        set(gca, 'XTick',datenum(1996:3:2023,0,0));
        set(gca, 'XTickLabel',yeara(:,3:4))
        axis([datenum(1995,12,1),datenum(2023,1,1),mineta(j), maxeta(j)])
        %axis([datenum(1995,12,1),datenum(2023,1,1),min(eta(:)), max(eta(:))])
        set(gca,'Box','on','LineWidth',4,'FontSize', 25)
        print('-depsc','-r70', [paperfd,'fig',pname{j},'rpts.eps'])

    end

end

if pricingperformance

    vca=100*[prctile(VC(:,:,1),[25,50,75])', prctile(VC(:,:,2),[25,50,75])'];
    disp('variance contribution');
    fprintf(1,' \n')
    for j=1:np
        fprintf(1, ['%20s '], pname{j});
        fprintf(1, [' &  %8.2f  &  %8.2f &  %8.2f &&  %8.2f  &  %8.2f &  %8.2f \\\\ \n'], vca(j,:)');
    end
    fprintf(1, ['%20s '], '$R^2$');
    fprintf(1, [' &  %8.2f  &  %8.2f &  %8.2f &&  %8.2f  &  %8.2f &  %8.2f \\\\ \n'], vca(np+1,:)');
    
    %pricing error in implied vol space
    eiv=reshape(evv,Tw,nm,nk,nc);
    ea1=reshape(nanmean(abs(eiv(:,:,:,1)))*1e2,nm,nk);
    figure(1);clf
    mesh(xm,ym,ea1,'LineWidth',4);
    xlabel('Put delta, %')
    ylabel('Maturity, Months')
    zlabel('Mean absolute error, %')
    axis([0,100,0,12,0,1])
    set(gca,'Box','on','LineWidth',4,'FontSize', 25)
    print('-depsc','-r70', [paperfd,'figMAE-JPY.eps'])
    
    ea2=reshape(nanmean(abs(eiv(:,:,:,2)))*1e2,nm,nk);
    figure(2);clf;
    mesh(xm,ym,ea2,'LineWidth',4);
    xlabel('Put delta, %')
    ylabel('Maturity, Months')
    zlabel('Mean absolute error, %')
    axis([0,100,0,12,0,1])
    set(gca,'Box','on','LineWidth',4,'FontSize', 25)
    print('-depsc','-r70', [paperfd,'figMAE-GBP.eps'])
    
    disp('mean absolute error');mae0=100*[nanmean(ea1(:)), nanmean(ea2(:))]


end

if returnfactormodel
        x=[FactorModel(:,:,1),FactorModel(:,:,2)];
        Tn=sum(isfinite(x));
        a=[52*nanmean(x),
            sqrt(52)*nanstd(x),
            sqrt(52)*nanmean(x)./nanstd(x),
            sqrt(Tn).*nanmean(x)./neweystd(x),
            ];
        
        a2=[100*nanmean(x),
            ];
        
        x=[RVC(:,:,1),RVC(:,:,2)];
        a3=100*nanmean(x);
        a4=[sum(a3(1:np+1)), sum(a3(np+2:end))];
        
        disp('Factor return behaviors');
        fprintf(1,' \n')
        for j=1:np+1
            fprintf(1, ['%20s '], pname{j});
            fprintf(1, [' &  %8.2f  &  %8.2f  &  %8.2f &  %8.2f &  %8.2f     &&  %8.2f  &  %8.2f  &  %8.2f  &  %8.2f  &  %8.2f    \\\\ \n'],    [a(:,j)',a3(j), a(:,np+2+j)',a3(np+1+j)]')
        end
        j=np+2;   fprintf(1,'\\\\ \n')
        fprintf(1, ['%20s '], '$R^2$');
        fprintf(1, [' &  ---  &  ---  &  --- &  --- & %8.2f    &&  ---  &  ---   &  --- &  --- & %8.2f     \\\\ \n'],    [a4']')
        
        
    fprintf(1,' \n');disp('factor return predictive regression');
    bbn2=NaN(np,6,nc);
    for c=1:nc
        for k=1:np
            y=52*FactorModel(:,k,c,1);
            x=RRP(:,k,c);
            ii=(isfinite(x+y));Tn=sum(ii);
            Xk=[ones(Tn,1), x(ii)];yk=y(ii);
            Bk=Xk\yk; ek=yk-Xk*Bk; R2k=1-mean(ek.^2)/var(yk);
            sbk=neweystd(ek)*sqrt(diag(inv(Xk'*Xk)));
            tbk=Bk./sbk;
            tbk2=(Bk-[0;1])./sbk;
            bbn2(k,:,c)=[Bk(1), tbk(1) Bk(2), tbk(2),tbk2(2), 100*R2k];
        end
    end
    fprintf(1,' \n')
    for j=1:np
        fprintf(1, [' %20s '], pname{j});
        fprintf(1, [' &  %8.2f   & ( %8.2f ) &&  %8.2f  & ( %8.2f  &  %8.2f ) &&  %8.2f  &&  %8.2f  & ( %8.2f )  &&  %8.2f  & ( %8.2f  &  %8.2f ) &&  %8.2f   \\\\ \n'], [bbn2(j,:,1),bbn2(j,:,2)]');
    end
    
    fprintf(1,' \n');disp('Investment performance');
    aan=NaN(np+1,6);
    aas=NaN(np+1,6);
    for c=1:2
        etav=NRP(:,:,c); %market pricing per notional
        rv=FPL(:,1:np,2,c); %future portfolio return per notional
  
        %cross-sectional normalization-- this is better as it does not
        %need to look forward, and it is easier to compare with the fixed
        %notional investment
        wv=etav./nanmean(abs(etav),2);
        trv=rv.*wv;
        ii=isfinite(sum(trv,2));
        Tn=sum(ii);
        x=[trv(ii,:), mean([trv(ii,:)],2)];
        aan(:,(c-1)*3+[1:3])=  [(52)*nanmean(x);   sqrt(52)*nanstd(x); sqrt(52)*nanmean(x)./nanstd(x)]';
        
        %static holding
        mr=-nanmean(reshape(PLD(:,:,:,c),Tw,nmk),2);
        x=[rv(ii,:),  mean(mr(ii,:),2)];
        aas(:,(c-1)*3+[1:3])=  [(52)*nanmean(x);   sqrt(52)*nanstd(x); sqrt(52)*nanmean(x)./nanstd(x)]';
        
        
    end
    
       fprintf(1,'fixed notional short \n')
    for j=1:np
        fprintf(1, [' %20s '], pname{j});
        fprintf(1, [' &  %8.2f  &  %8.2f  &  %8.2f    &&  %8.2f  &  %8.2f  &  %8.2f     \\\\ \n'], [aas(j,:)]');
    end
    fprintf(1,' \n');j=np+1;
        fprintf(1, [' %20s '], '\\Market portfolio');
        fprintf(1, [' &  %8.2f  &  %8.2f  &  %8.2f    &&  %8.2f  &  %8.2f  &  %8.2f     \\\\ \n'], [aas(j,:)]');

        fprintf(1,' market timing\n')
    for j=1:np
        fprintf(1, [' %40s '], pname{j});
        fprintf(1, [' &  %8.2f  &  %8.2f  &  %8.2f    &&  %8.2f  &  %8.2f  &  %8.2f     \\\\ \n'], aan(j,:)');
    end
    fprintf(1,' \n');
    j=np+1;
    fprintf(1, [' %40s '], '\\Risk premium harvesting');
    fprintf(1, [' &  %8.2f  &  %8.2f  &  %8.2f    &&  %8.2f  &  %8.2f  &  %8.2f     \\\\ \n'], aan(j,:)');

end




return




if summarystats
    x1=reshape(IV(:,:,:,1), T,nmk)*100;
    x2=reshape(IV(:,:,:,2), T,nmk)*100;
    figure(1)
    clf
    plot(ddate,x1,'LineWidth',2)
    xlabel('Year');
    ylabel('Implied volatility, %');
    set(gca, 'YTick',[0:20:100]);
    set(gca, 'XTick',datenum(1996:3:2023,0,0));
    set(gca, 'XTickLabel',yeara(:,3:4))
    axis([datenum(1995,12,1),datenum(2023,1,1),0,60])
    % grid
    set(gca,'Box','on','LineWidth',4,'FontSize', 25)
    print('-depsc','-r70', [paperfd,'figIVats_jpy.eps'])
    
    figure(2)
    clf
    plot(ddate,x2,'LineWidth',2)
    xlabel('Year');    ylabel('Implied volatility, %');
    set(gca, 'YTick',[0:20:100]);
    set(gca, 'XTick',datenum(1996:3:2023,0,0));
    set(gca, 'XTickLabel',yeara(:,3:4))
    axis([datenum(1995,12,1),datenum(2023,1,1),0,60])
    %grid
    set(gca,'Box','on','LineWidth',4,'FontSize', 25)
    print('-depsc','-r70', [paperfd,'figIVats_gbp.eps'])
    
    %term structure
    ts1=100*(IV(:,5,3,1)-IV(:,1,3,1));
    ts2=100*(IV(:,5,3,2)-IV(:,1,3,2));
    
    sk1=100*(IV(:,2,5,1)-IV(:,2,1,1));
    sk2=100*(IV(:,2,5,2)-IV(:,2,1,2));
    
    figure(3)
    clf
    plot(ddate,ts1,ddate,ts2,'-.',ddate,ddate*0,':','LineWidth',4)
    xlabel('Year');
    ylabel('Implied volatility term structure, %');
    set(gca, 'YTick',[-20:5:5]);
    set(gca, 'XTick',datenum(1996:3:2023,0,0));
    set(gca, 'XTickLabel',yeara(:,3:4))
    axis([datenum(1995,12,1),datenum(2023,1,1),-20,5])
    %grid
    set(gca,'Box','on','LineWidth',4,'FontSize', 25)
    print('-depsc','-r70', [paperfd,'figIVtsts.eps'])
    
    figure(4)
    clf
    plot(ddate,sk1,ddate,sk2,'-.',ddate,ddate*0,':','LineWidth',4)
    xlabel('Year');
    ylabel('Implied volatility skew, %');
    set(gca, 'YTick',[-20:5:20]);
    set(gca, 'XTick',datenum(1996:3:2023,0,0));
    set(gca, 'XTickLabel',yeara(:,3:4))
    axis([datenum(1995,12,1),datenum(2023,1,1),-20,20])
    %grid
    set(gca,'Box','on','LineWidth',4,'FontSize', 25)
    print('-depsc','-r70', [paperfd,'figIVskts.eps'])
    
    
    
    x=nanmean(CSMoments); y=reshape(nanmean(CSMoments(:,:,:,4,:),2),Tw,4,2);
    a3=reshape(x(1,:,4,4,1),nm,nk);
    figure(5);clf;
    plot(matm, a3(:,1),'-',matm, a3(:,2),'--',matm, a3(:,3),'-.',matm, a3(:,4),'o:',matm, a3(:,5),':','LineWidth',4,'MarkerSize',20)
    xlabel('Maturity, Months')
    ylabel('\omega')
    %axis([0,100,0,12,-50,50])
    set(gca,'Box','on','LineWidth',4,'FontSize', 25)
    print('-depsc','-r70', [paperfd,'figcsomega-JPY.eps'])
    
    a3=reshape(x(1,:,4,4,2),nm,nk);
    figure(6);clf;
    plot(matm, a3(:,1),'-',matm, a3(:,2),'--',matm, a3(:,3),'-.',matm, a3(:,4),'o:',matm, a3(:,5),':','LineWidth',4,'MarkerSize',20)
    xlabel('Maturity, Months')
    ylabel('\omega')
    %axis([0,100,0,12,-50,50])
    set(gca,'Box','on','LineWidth',4,'FontSize', 25)
    print('-depsc','-r70', [paperfd,'figcsomega-GBP.eps'])
    
    x=CSMoments(:,11:15,4,4,1)'; %almost always downward sloping
     figure(7);clf;
    plot(matm, log(x),'LineWidth',2,'MarkerSize',20)
    xlabel('Maturity, Months')
    ylabel('\omega')
    %axis([0,100,0,12,-50,50])
    set(gca,'Box','on','LineWidth',4,'FontSize', 25)
    print('-depsc','-r70', [paperfd,'figcsomegatsts-JPY.eps'])

    
    
    % mom(:,:,h)=[hmt,hvt,hgt,hot];

    y=reshape(CSMoments(:,13,:,4,:),Tw,4,2);v=reshape(y(:,2,:),Tw,2);    o=reshape(y(:,4,:),Tw,2);


    figure(7)
    clf
    plot(wdate,y(:,4,1),wdate,y(:,4,2),'-.','LineWidth',4)
    xlabel('Year');
    ylabel('\omega_t');
    %set(gca, 'YTick',[-20:5:20]);
    set(gca, 'XTick',datenum(1996:3:2023,0,0));
    set(gca, 'XTickLabel',yeara(:,3:4))
    axis([datenum(1995,12,1),datenum(2023,1,1),0,3])
    %grid
    set(gca,'Box','on','LineWidth',4,'FontSize', 25)
    print('-depsc','-r70', [paperfd,'figcsomegats.eps'])
    
    figure(8)
    clf
    plot(wdate,y(:,2,1),wdate,y(:,2,2),'-.','LineWidth',4)
    xlabel('Year');
    ylabel('v_t');
    %set(gca, 'YTick',[-20:5:20]);
    set(gca, 'XTick',datenum(1996:3:2023,0,0));
    set(gca, 'XTickLabel',yeara(:,3:4))
    axis([datenum(1995,12,1),datenum(2023,1,1),0,0.12])
    %grid
    set(gca,'Box','on','LineWidth',4,'FontSize', 25)
    print('-depsc','-r70', [paperfd,'figcsvts.eps'])
    
    figure(9)
    clf
    plot(wdate,BO(:,1),wdate,BO(:,2),'-.',wdate, wdate*0,':','LineWidth',4)
    xlabel('Year');
    ylabel('\rho_t');
    %set(gca, 'YTick',[-20:5:20]);
    set(gca, 'XTick',datenum(1996:3:2023,0,0));
    set(gca, 'XTickLabel',yeara(:,3:4))
    axis([datenum(1995,12,1),datenum(2023,1,1),-1,1])
    %grid
    set(gca,'Box','on','LineWidth',4,'FontSize', 25)
    print('-depsc','-r70', [paperfd,'figcsrhots.eps'])
    
    
    
    
    mPL=-[reshape(nanmean(PLD(:,:,:,1))*52,nm,nk),reshape(nanmean(PLD(:,:,:,2))*52,nm,nk)];
    sPL=[reshape(nanstd(PLD(:,:,:,1))*sqrt(52),nm,nk),reshape(nanstd(PLD(:,:,:,2))*sqrt(52),nm,nk)];
    iPL=mPL./sPL;
    
    disp('Mean delta-hedged return, annaulized %')
    fprintf(1,' \n')
    fprintf(1, ['%8.0f  ', repmat(' &  %8.2f  ',1,5), ' & ',repmat(' &  %8.2f  ',1,5),' \\\\ \n'], [matm,mPL]');
    disp('Vol')
    fprintf(1,' \n')
    fprintf(1, ['%8.0f  ', repmat(' &  %8.2f  ',1,5), ' & ',repmat(' &  %8.2f  ',1,5),' \\\\ \n'], [matm,sPL]');
    disp('IR')
    fprintf(1,' \n')
    fprintf(1, ['%8.0f  ', repmat(' &  %8.2f  ',1,5), ' & ',repmat(' &  %8.2f  ',1,5),' \\\\ \n'], [matm,iPL]');
    disp('equal average')
    x=[nanmean(reshape(PLD(:,:,:,1),Tw,nmk),2),nanmean(reshape(PLD(:,:,:,2),Tw,nmk),2)];
    [-nanmean(x)*52; nanstd(x)*sqrt(52); -nanmean(x)*sqrt(52)./nanstd(x)]
    
    x(~isfinite(x))=0;cpl=cumsum(-x);
    figure(9);
    clf
    plot(wdate,cpl(:,1),wdate,cpl(:,2),'-.','LineWidth',4)
    xlabel('Year');    ylabel('Cumulative return, %');
    set(gca, 'YTick',[-10:2:20]);
    set(gca, 'XTick',datenum(1996:3:2023,0,0));
    set(gca, 'XTickLabel',yeara(:,3:4))
    axis([datenum(1995,12,1),datenum(2023,1,1),-2,12])
    %grid
    set(gca,'Box','on','LineWidth',4,'FontSize', 25)
    print('-depsc','-r70', [paperfd,'figEWCPLts.eps'])
    
end





return
figure(27)
clf;
plot(dd,100*AR2(:,1), dd,100*AR2(:,2),'-.','LineWidth',4)
xlabel('Year','FontSize',25);
ylabel('Adjusted R^2, %','FontSize',25);
set(gca, 'YTick',[0:20:100]);
set(gca, 'XTick',datenum(1996:3:2023,0,0));
set(gca, 'XTickLabel',yeara(:,3:4))
axis([datenum(1995,12,1),datenum(2023,1,1),0,110])
set(gca,'Box','on','LineWidth',3,'FontSize', 25)
print('-depsc','-r70', [paperfd,'figAR2ts.eps'])

figure(28)
clf;
plot(dd,100*R22(:,1), dd,100*R22(:,2),'-.','LineWidth',4)
xlabel('Year','FontSize',25);
ylabel('R^2, %','FontSize',25);
set(gca, 'YTick',[0:20:100]);
set(gca, 'XTick',datenum(1996:3:2023,0,0));
set(gca, 'XTickLabel',yeara(:,3:4))
axis([datenum(1995,12,1),datenum(2023,1,1),0,110])
set(gca,'Box','on','LineWidth',3,'FontSize', 25)
print('-depsc','-r70', [paperfd,'figR2ts.eps'])

x=-[FPL(:,np+1,1,1),FPL(:,np+1,1,2),FPL(:,np+1,2,1),FPL(:,np+1,2,2)];
ii=isfinite(sum(x,2));
cx=cumsum(x(ii,:));ddi=dd(ii);

figure(29)
clf;
plot(ddi,cx(:,1),ddi,cx(:,2),'-.','LineWidth',4)
xlabel('Year','FontSize',25);
ylabel('Cumulative P&L','FontSize',25);
legend('JPY','GBP','Location','NorthWest')
set(gca, 'YTick',[0:40:160]);
set(gca, 'XTick',datenum(1997:3:2023,0,0));
set(gca, 'XTickLabel',yeara(:,3:4))
axis([datenum(1996,12,1),datenum(2023,1,1),0,160])
set(gca,'Box','on','LineWidth',3,'FontSize', 25)
print('-depsc','-r70', [paperfd,'figcPL6ts.eps'])



a1=reshape(nanmean((Iw(:,:,:,1))),nm,nk)*100;
figure(1);clf
mesh(xm,ym,a1,'LineWidth',4);
xlabel('Put delta, %')
ylabel('Maturity, Months')
zlabel('Implied volatility, %')
axis([0,100,0,12,10,14])
set(gca,'Box','on','LineWidth',4,'FontSize', 25)
print('-depsc','-r70', [paperfd,'figIV-JPY.eps'])


a2=reshape(nanmean((Iw(:,:,:,2)))*1e2,nm,nk);
figure(2);clf;
mesh(xm,ym,a2,'LineWidth',4);
xlabel('Put delta, %')
ylabel('Maturity, Months')
zlabel('Implied volatility, %')
axis([0,100,0,12,8,12])
set(gca,'Box','on','LineWidth',4,'FontSize', 25)
print('-depsc','-r70', [paperfd,'figIV-GBP.eps'])





figure(7)
clf;
plot(dd,y(:,1,1), dd,y(:,1,2),'-.',dd,dd*0,':','LineWidth',4)
xlabel('Year','FontSize',25);
ylabel('\mu_t','FontSize',25);
set(gca, 'YTick',[-4:2:6]);
set(gca, 'XTick',datenum(1996:3:2023,0,0));
set(gca, 'XTickLabel',yeara(:,3:4))
axis([datenum(1995,12,1),datenum(2023,1,1),-4,6])
set(gca,'Box','on','LineWidth',3,'FontSize', 25)
print('-depsc','-r70', [paperfd,'figcsmuts.eps'])


figure(8)
clf;
plot(dd,y(:,4,1), dd,y(:,4,2),'-.','LineWidth',4)
xlabel('Year','FontSize',25);
ylabel('\omega_t','FontSize',25);
set(gca, 'YTick',[0:.5:2]);
set(gca, 'XTick',datenum(1996:3:2023,0,0));
set(gca, 'XTickLabel',yeara(:,3:4))
axis([datenum(1995,12,1),datenum(2023,1,1),0,2])
set(gca,'Box','on','LineWidth',3,'FontSize', 25)
print('-depsc','-r70', [paperfd,'figcsomegats.eps'])

figure(9)
clf;
plot(dd,y(:,2,1), dd,y(:,2,2),'-.','LineWidth',4)
xlabel('Year','FontSize',25);
ylabel('v_t','FontSize',25);
set(gca, 'YTick',[0:.02:.1]);
set(gca, 'XTick',datenum(1996:3:2023,0,0));
set(gca, 'XTickLabel',yeara(:,3:4))
axis([datenum(1995,12,1),datenum(2023,1,1),0,.1])
set(gca,'Box','on','LineWidth',3,'FontSize', 25)
print('-depsc','-r70', [paperfd,'figcsvts.eps'])

figure(10);clf;
plot(dd, BO(:,1),dd,BO(:,2),'-.',dd,dd*0,':','LineWidth',4)
xlabel('Year','FontSize',25);
ylabel('\rho_t','FontSize',25);
set(gca, 'YTick',[-1:.5:1]);
set(gca, 'XTick',datenum(1996:3:2023,0,0));
set(gca, 'XTickLabel',yeara(:,3:4))
axis([datenum(1995,12,1),datenum(2023,1,1),-.5,.5])
set(gca,'Box','on','LineWidth',3,'FontSize', 25)
print('-depsc','-r70', [paperfd,'figrhots.eps'])




tt=[1:T]; pv=[10:90]';np=length(pv);
xh=NaN(np,2,np);yh=xh;
for j=1:np
    eta=[RRP(:,j,1), RRP(:,j,2)];
    figure(12+j);clf;
    plot(dd, eta(:,1),dd,eta(:,2),'-.',dd,dd*0,':','LineWidth',4)
    xlabel('Year','FontSize',20);
    ylabel(xpname{j},'FontSize',20);
    set(gca, 'XTick',datenum(1996:3:2023,0,0));
    set(gca, 'XTickLabel',yeara(:,3:4))
    axis([datenum(1995,12,1),datenum(2023,1,1),min(eta(:)), max(eta(:))])
    set(gca,'Box','on','LineWidth',4,'FontSize', 25)
    print('-depsc','-r70', [paperfd,'fig',pname{j},'rpts.eps'])
    
    raw=0;
    if raw
        eta=[RRP(:,j,1), RRP(:,j,2)];
        pr=-[FPL(:,j,1,1),FPL(:,j,1,2)]; %risk premium raw
    else
        eta=[RP(:,j,1), RP(:,j,2)];
        pr=-[FPL(:,j,2,1),FPL(:,j,2,2)]; %risk premium per risk
    end
    
    for c=1:2
        xh(:,c,j)=prctile(eta(:,c),pv);
        yh(:,c,j)=llregressr(eta(:,c),pr(:,c),xh(:,c,j),4)*252;
    end
    
end
for j=1:np
    figure(17+j)
    clf;
    plot(xh(:,1,j), yh(:,1,j), xh(:,2,j), yh(:,2,j),'-.',0,0,'o','MarkerSize',35,'LineWidth',5)
    xlabel(xpname{j},'FontSize',20);
    ylabel('Portfolio return, %','FontSize',20);
    grid
    set(gca,'Box','on','LineWidth',4,'FontSize', 25)
    print('-depsc','-r70', [paperfd,'fig',pname{j},'PLfit.eps'])
end

figure(25)
clf;
plot(xh(:,1,1), yh(:,1,1), xh(:,2,1), yh(:,2,1),'-.',0,0,'o', ...
    xh(:,1,2), yh(:,1,2), xh(:,2,2), yh(:,2,2),'-.', ...
    'MarkerSize',35,'LineWidth',5)
xlabel(['\eta_V, \eta_\mu'],'FontSize',20);
ylabel('Portfolio return, %','FontSize',20);
grid
set(gca,'Box','on','LineWidth',4,'FontSize', 25)
print('-depsc','-r70', [paperfd,'figVegaTrendPLfit.eps'])

